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ABSTRACT 

Mergers of two carbon-oxygen white dwarfs have long been suspected to be progen- 
itors of Type la Supernovae. Here we present our modifications to the cosmological 
smoothed particle hydrodynamics code Gadget to apply it to stellar physics includ- 
ing but not limited to mergers of white dwarfs. We demonstrate a new method to map 
a one-dimensional profile of an object in hydrostatic equilibrium to a stable particle 
distribution. We use the code to study the effect of initial conditions and resolution on 
the properties of the merger of two white dwarfs. We compare mergers with approxi- 
mate and exact binary initial conditions and find that exact binary initial conditions 
lead to a much more stable binary system but there is no difference in the properties of 
the actual merger. In contrast, we find that resolution is a critical issue for simulations 
of white dwarf mergers. Carbon burning hotspots which may lead to a detonation in 
the so-called violent merger scenario emerge only in simulations with sufficient res- 
olution but independent of the type of binary initial conditions. We conclude that 
simulations of white dwarf mergers which attempt to investigate their potential for 
Type la supernovae should be carried out with at least 10 6 particles. 

Key words: stars: supernovae: general - hydrodynamics - binaries: close - methods: 
numerical 



1 INTRODUCTION 

Mergers of two white dwarfs have first been proposed as 
progenitor systems of Type la Supernovae (SNe la) by |Iben| 
"fe Tutukov| ( |1984| ) and |Webbink| ( |1984| ). More than 25 years 
later, however, we are neither able to confirm them nor rule 
them out theoretically nor observationally. 

On the theory side there is a long tradition of numerical 
simulations of white dwarfs mergers to determine the fate of 
those systems after the merger. With time, these simulations 
have become better and more and more sophisticated both 
in terms of the numerical resolution used and the treatment 
of the input physics. 

The merger of two white dwarfs is an inherently three- 
dimensional problem. There is no intrinsic symmetry that 
can be exploited to simulate it in less then three dimensions. 
Starting with the pioneering work by Benz et al. ([1990 ) most 
simulations of white dwarf mergers used smoothed particle 
hydrodynamics (SPH) codes. This results from several ad- 
vantages SPH codes have compared to finite volume codes 
for this particular problem, including the fact that SPH con- 
serves angular momentum very well, does not need a descrip- 



tion of the volume surrounding an object and offers a simple 
way of implementing additional physics. Nevertheless, there 



are a few attempts to use finite volume codes (D'Souza et al. 



2006 



Motl et al.|2007| . A problem of these simulations, how- 



ever, has been shown to be the use of a polytropic equation 
of state which suppresses shocks and leads to an artificial 
increase of the orbital separation ( Dan et al.||2011 ). Quali- 
tatively, their results were not so different from |Benz et al.] 
(1990), although they used a resolution of only 7000 parti- 
cles. They showed that the merger of two white dwarfs will 
lead to the destruction of the secondary, less massive, white 
dwarf. Its material forms a hot envelope and an accretion 
disk around the remaining primary white dwarf. 



The result of the long-term evolution of such a merger 
remnant is then likely to be the conversion of the central core 
into an oxygen-neon white dwarf by a slowly inward prop- 



agating flame (e.g. |Nomoto &; Iben 1985 Saio & Nomoto 



1998). If this oxygen- neon white dwarf accretes enough ma- 



terial from its envelope and the accretion disk to reach the 
Chandrasekhar mass, it will undergo an accretion-induced 
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collapse, rather than a thermonuclear explosion (Nomoto & 
Kondo|[l99T ) 



Following iBenz et all dl990L merger simulations have 
been r epeat ed with more and more particles (i.e . |Guerrero | 
et al.| \2QQ4\ using 4 x 10 * particles, |Yoon et al | fl2007| ) us- 
ing 2 x 10 particles, and Loren-Aguilar et al. ( 2Q09| ) using 
4 x 10 5 particles) confirming previous results. Starting with 
Guerrero et al.| ( |2004| ) a nuclear reaction network has been 
coupled to the equations of hydrodynamics to account for 
the energy release from nuclear reactions during the merger. 

Only recently, indications have been found that those 
results may not be the full picture. Pakmor et al. (2010) who 
simulated the merger of two equal-mass white dwarfs with a 
mass ratio of one with an unprecedented resolution of 2 x 10 6 
particles found that in the final phase of the merger, just 
before the disruption of the secondary white dwarf, a hotspot 
develops on the surface of the primary white dwarf. In this 
hotspot carbon burning starts and a detonation may form 
that directly consumes the whole merging binary system 
leading to a thermonuclear explosion. 

This violent merger scenario, of course, stands and falls 
with the assumption that a detonation forms. Further work 
showed that these hotspots form even with lower mass ra- 



tios of the merging white dwarfs, down to 0.8 (Pakmor et al 



2011 ). Doubts on the formation of those hotspots have been 
raised by Dan et al. ( 2011| who introduced a method to ob- 
tain more relaxed ("exact") binary initial conditions. They 
found that this leads to a significantly more stable binary 
system and without hotspots at a resolution of 2 x 10 5 par- 
ticles. In contrast, |Raskin et al.| (2012) and Pakmor et al. 



( |2012| found hotspots to be present in merger simulations 
with exact binary initial conditions and a resolution of 10 6 
particles and more. Note, however, that all these simulations 
have been run with different codes, different parameters, and 
at least partially different input physics. 

In this paper, we therefore describe in detail our 
methodology, which is based on the Gadget code. While 
we emphasize that the modifications of the code, that orig- 
inally was designed to address problems in cosmology, are 
useful for a potentially wide range of problems in stellar as- 
trophysics, we use our implementation to test the role of 
initial conditions and resolution in mergers of white dwarfs 
here. 

In Sec.[2]we describe in detail the modification we made 
to the Gadget code and test our new timestepping mecha- 
nism in Sec. [3] We present our new nuclear reaction network 
in Sec. We discuss our new method to create initial con- 
ditions for stars in SPH in Sec. \S\ and demonstrate that it 
works well. In addition, we describe our implementation to 
create "exact" binary initial conditions following |Dan et al.| 
( 2011 ). We then show results of simulations of merging white 



dwarfs in Sec. [5] studying in particular the influence of bi- 
nary initial conditions and resolution on the merger of two 
massive white dwarfs using the same code and parameters 
for all simulations. We finish with a summary of the paper 
and some conclusions in Sec. 



2 CODE MODIFICATIONS 

A detailed description of the Gadget code can be found 



include additional physics as well as technical changes to 
adapt the code to different conditions. The changes are as 
follows. 

(i) Gravitational softening 
To calculate the gravitational forces we use the standard 
tree algorithm implemented in Gadget. Instead of a 
fixed gravitational softening length as usually applied in 
cosmological simulations, we use the individual smoothing 
lengths of gas particles as their gravitational softening 



(Bate & Burkert 1997). This choice significantly improves 
the stability of objects in hydrostatic equilibrium. How- 
ever, it can lead to errors in the total energy budget of 
the simulation when the smoothing length of a particle 
changes, as the gravitational softening and therefore the 
local gravitational potential also changes. As described by 



Price & Monaghan (2007), it is possible to compensate 



for this by adding extra terms to the evaluation of the 
gravitational force (for an implementation into the Gadget 
code see |Iannuzzi fc Dolag (2011)). The improved energy 
conservation, however, comes at the cost of additional 
noise in the velocity field. Therefore, we refrain from using 
it. For most applications the errors introduced are small 
anyway, because the smoothing lengths of the particles do 
not change significantly during the simulation. 

(ii) Wakeup mechanism 
In the Gadget code, individual time steps are assigned to all 
particles, which depend only on the local conditions near the 
particles. The particles are then evolved on their own time 
steps. While this greatly reduces the computational costs to 
run a simulation, it can cause severe problems if very fast 
particles run into a medium in which the sound speed is 
much smaller then the velocity of the fast particles. In this 
case, a fast particle on a small time step can pass through a 
particle on a much larger time step without being noticed. 
As this behavior can obviously lead to completely unphys- 
ical results, we implement a "wakeup mechanism" for the 
time-stepping. Similar to the method proposed by Sait oh &\ 
Makino ( 2009 ) , it should activate inactive particles as other 



particles approach them which evolve on much shorter time 
steps. 

The hydro dynamical time step of an SPH particle is calcu- 
lated as 

Ccourant hi 



AU 



signal \ 



(1) 



where C CO urant is the Courant factor, hi is the smoothing 
length of the particle and the denominator finds the maxi- 
mum of the signal velocities v^ gnal from particle i and all its 



neighbors j as defined in equation 13 of Springel (2005). The 



in Springel (2005) and Springel et al. ( 2001| ). Our changes 



resulting maximum signal velocity of a particle is stored. In 
each time step we check for all inactive neighbors of all ac- 
tive particles whether their stored maximum signal velocity 
is smaller by some factor than the new signal velocity cal- 
culated for the active particle and the inactive neighbor, 

(2) 

If fulfilled, we flag the particle to be woken up. Our usual 
choice for the wakeup factor w; is 4.1; thus an active particle 
can be in range of an inactive particle for three time steps 
at most. The check can be done very efficiently, since the 
active particle has to loop over all its neighbors anyway to 



maxsignal ^ signal 

v i s < w ■ . 
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calculate the local pressure force. After finishing the time 
step, we change the time step of all particles that have been 
flagged to wake up such that they become active in the next 
possible time step. When a particle was active in a time 
step, its properties had been predicted half a time step into 
the future using its actual rates of change of these proper- 
ties. Thus, this extrapolation needs to be corrected when 
we shorten the time step. The correction for the difference 
between the estimated time step and the time step the par- 
ticle really experienced is done consistently for all quantities. 

(hi) Energy equation 
In contrast to the original implementation of Gadget that 
uses the entropy equation we solve the energy equation. 
This is convenient choice because evolving the internal en- 
ergy simplifies combining hydrodynamics with the nuclear 
reaction network. 

(iv) Equation of state 

The equation of state (EoS) captures all intrinsic properties 
of the material. It is used to calculate local pressure 
and speed of sound from the primary thermodynamical 
variables evolved in the code. The original Gadget code 
only implements an ideal gas equation of state. We replace 
it with the Helmholtz-EoS flTimmes fc Swesty| |2000| ). 
This EoS describes an arbitrarily degenerate, arbitrarily 
relativistic electron-positron gas together with an ideal gas 
of completely ionized ions. It also includes radiation from a 
black body with the local gas temperature. Since internal 
energy is our thermodynamic variable of choice, most calls 
of the EoS have the internal energy as input. Because it is 
in general not possible to invert the EoS in a closed form, 
we need to iterate on the temperature first, using a Newton 
method. Afterwards all other quantities are calculated from 
temperature, density and composition. Due to numerical 
errors (e.g. when the value of the internal energy becomes 
smaller than the minimum energy of the degenerate elec- 
tron), it may not be possible to find a valid temperature 
for a given internal energy. In this rare case, or when the 
temperature drops below the minimum temperature tabu- 
lated by the EoS, we assign a temperature of 1000K to the 
particle. All other thermodynamical quantities apart from 
the internal energy are then calculated for this temperature. 

(v) Nuclear reaction network 

Because SPH is a purely Lagrangian method, the composi- 
tion of a particle can only change by nuclear reactions. The 
nuclear network is coupled to hydrodynamics in an opera- 
tor splitting approach. After a hydrodynamics timestep, the 
nuclear network is integrated for all active particles for the 
duration of their timestep. During this integration, the den- 
sity is assumed to be constant. Changes of the temperature 
as a result of energy release or consumption by nuclear re- 
actions are taken into account in the network. The nuclear 
reaction network and its integration are described in detail 
in Section [3] We restrict evolving the nuclear network to par- 
ticles with temperatures higher than 10 7 K. For most of our 
applications this is a very conservative estimate, as carbon 
burning starts only at around 10 9 K. 

After evolving the nuclear network for the duration of the 
hydrodynamical time step, the composition of the particle 
is changed according to the nuclear reaction network. From 
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Figure 1. Radial density profiles of two simulation runs of the 
Sedov problem with global timesteps and local timesteps with the 
wakeup mechanism, respectively. The analytical solution is shown 
for comparison. 



the rate of change of the abundances of the different species 
we then calculate the amount of energy that is released or 
consumed by nuclear reactions. 



n a M 3 c 2 



dt 



(3) 



Here, Na is Avogadro's constant, MjC 2 is the rest mass en- 
ergy of species j, and Yj is its number fraction. This change 
of energy is included as a source term in the energy equation. 



3 WAKEUP TEST: SEDOV PROBLEM WITH 
LOCAL TIMESTEPPING 

To check the accuracy of our wakeup mechanism we apply 
it to the Sedov problem (Sedov 1959}, which is a simple 



point explosion leading to a strong shock wave expanding 
in a surrounding medium. We start from glass initial condi- 
tions flWhite||1994 ) in a periodic box of size 10 in arbitrary 
units that contains 10 6 particles of with a mass of 10~ 3 each. 
Then we inject an internal energy of 10 5 equally into the 33 
particles closest to the center of the box. The internal en- 
ergy of the surrounding particles is chosen to be 10 -5 . For 
this configuration there is a large spread of local timesteps. 
Particles which are influenced by one of the hot particles 
(i.e. one of them is closer than their smoothing length) are 
evolved on a timestep of 10 -5 . All other particles, however, 
start with a timestep of 5 x 10 2 , more than a million times 
larger. Therefore, these particles obviously decouple com- 
pletely from the hot particles and their close neighbours. 
Although the conditions we use are quite extreme, the same 
problem also occurs for real problems, i.e. the interaction of 
the eject a of a supernova explosion with a companion star 
in a binary system (e.g. |Pakmor et al. 2008). One way to 
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circumvent this problem is to use the same global timestep 
for all particles, which guarantees correct behavior. How- 
ever, it makes the simulation significantly more expensive, 
because we then have to evolve all particles on a very small 
timestep from the beginning rather than from the time when 
the shock actually hits them. Another approach is discussed 
in Sec. [2] which reduces the timestep of a particle as soon as 
the shock touches it. Fig. [I] compares the radial density pro- 
file for the Sedov problem for a global timestep and a local 
timestep combined with this wakeup mechanism. It clearly 
shows that using local timesteps plus the wakeup mecha- 
nism recovers exactly the same result we get from using a 
global timestep. 



4 THE NUCLEAR NETWORK 

The nuclear reaction network calculates the change of the 
composition with time. The change of the abundance Y$ of 
one species depends on the reaction rates IZk that destroy 
or create nuclei of this species. 



Y i = J2^k (p,T,Y) 



(4) 



Each reaction rate depends on the local density, tempera- 
ture and abundance of the react ants. The reaction rates are 
tabulated and taken from the latest (2009) release of the 
REACLIB database ( |Rauscher & Thielemann||2000| ), which 
includes experimental data as well as theoretical rates when 
experimental data are not available. Additionally, the weak 
interaction rates of Langanke & Martmez-Pinedo (2001 ) are 
included. Neutrino losses are neglected. 

Under the conditions relevant for nucleosynthesis atoms 
are usually fully ionized. Nevertheless, the assumption of 
bare nuclei colliding with each other is not valid, because the 
electron background shields part of the charge of the nuclei. 
Thus, since the Coulomb repulsion is reduced, the nuclei 
can interact more easily and the reaction rate increases. We 
treat the effects of electron screening on the reaction rates 



as described by Wallace et al. (1982). This description dis- 



criminates between regimes of strong ( Alastuey &; Jancovici| 
|1978| and weak ( |Dewitt et al]|1973t screening and an in- 



termediate region. In contrast to the original paper, we use 
a different interpolation method in the intermediate regime 
and a corrected coefficient for strong screening, both taken 
from I Wallace et al.| p983) . 



Since the reaction rates depend on temperature to a 
high power and temperature can change quickly due to the 
energy release of nuclear reactions, the temperature has to 
be evolved along with the abundances for correct integra- 
tion. This is taken care of by an additional equation for the 
evolution of temperature 



T= E 



dT 
dE 



p=const 



dT 
Wi' 



(5) 



It should be noted that T is completely determined by the 
abundances Yi and the density, if the energy release is given 
by equation (|3| . Therefore it is not necessary to include it ex- 
plicitly as an additional variable. The problem with this ap- 
proach is that the Jacobian of the right side of equation Q 
would be a dense matrix. By separating the temperature 
evolution into an additional equation, the Jacobian becomes 



sparse, with non-zero entries at the respective input and 
output nuclei of the corresponding reactions. Therefore, we 
choose to include the extra equation for T to be able to use 
significantly faster algorithms for inverting the Jacobian. 

We assume the density to be constant over one timestep. 
Once this approximation is not valid anymore, an evolution 
equation for the density can be implemented the same way 
as for the temperature. Some kind of thermodynamic con- 
straint (e.g. constant pressure) is needed in this case. 

If electron screening is accounted for, the same method 
is used for the quantity ^^p ecies Zj 2 Yj on which the 
screening factors depend. The same would apply to Y e = 
y^jspecies ZjYj but because Y e only changes slowly over a 
typical timestep it can safely be assumed be constant for 
the computation of the Jacobian. 

Thus, we end up with a system of N + 1 (N + 2, with 
screening) ordinary differential equations for N + 1 (N + 2, 
with screening) variables in a nuclear network of N species. 
Because this system of equations is very stiff, it needs to 
be integrated implicitly. To this end we apply the variable- 



order Bader-Deuflhard method (Bader & Deuflhard 1983) as 



suggested by |Timmes| ( |1999| ) . Depending on the number of 
nuclei we use a full direct matrix solver (for small networks) 
or the sparse matrix solver SuperLU ( Demmel et al.|1 999). 



5 INITIAL CONDITIONS FOR COMPACT 
OBJECTS 

Models of stars and compact objects usually assume hydro- 
static equilibrium and spherical symmetry. When we want 
to simulate such objects in Gadget, we therefore have to 
start from one-dimensional profiles of density, composition, 
pressure, etc. in hydrostatic equilibrium. We then have to 
find a way to construct a stable three-dimensional particle 
distribution, that reproduces these profiles. Because the lo- 
cal density is coupled to the local particle distribution, noise 
in the density estimate of any non-trivial initial particle con- 
figuration cannot be avoided. It is essential to reduce this 
noise to maintain the initial configuration as accurate as 
possible. 



5.1 Setup 

The idea of our method to construct the particle distribution 
of a single star is to divide the star into spherical shells and 
these shells into subvolumes of roughly cubical size to which 
we attribute one SPH particle each. In addition, each of 
the subvolumes should contain the same mass, to fulfill the 
constraint of equally massive particles. With this method we 
obtain a rather regular particle distribution that resembles 
any given radial density profile with only very small noise 
compared to a random sampling of the density profile. 

In our implementation we use the healpix library 
(Gor ski et al.|2005 ) to tessellate the surface of a sphere into 
12 • N z approximately quadrat ical pieces of the same area 
where the index N is a non-negative integer. We then con- 
struct the star out of several shells starting from the center. 
The width of the shells is chosen such that it is of about 
the same size as the edges of the pieces. As we show below, 
this condition can be satisfied using a simple method to de- 
termine the width of the shells. We express the constraints 
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Figure 2. Comparison of the initial conditions of a I.OMq white dwarf generated by the healpix method (see text). Black points show 
the SPH density estimate of all particles depending on their radial coordinate. The red line shows the initial one-dimensional density 
profile. Left and right column show the same data, but use a linear/logarithmic scale for the density. The top row shows the initial setup, 
the bottom row shows the configuration after completion of the relaxation. 



using parameters m and 712, which are equal if both con- 
straints are fulfilled. Then N = round(ni) = round(ri2) is 
used to tesselate the shell. 

Let us assume we know the lower radius ri OW er of a shell 
and want to find its upper radius r upP er- For a given upper 
radius, the width of a shell and its mass are given by 



ttshell — Supper Slower 

and 

/'''"upper 

m s heii = 47T / p(r)r 2 dr, 

J n ower 



(6) 



(7) 



respectively. The mass increases with increasing r upP er. Fix- 
ing the uniform particle mass in the beginning, we can cal- 
culate the number of particles that need to be placed in this 
shell to 



^particles 



TTlshell 



'^particle 

This is equivalent to an index m of 



m — 



^shell 



12 • m 



particle 



(8) 



(9) 



which increases with increasing width of the shell. 

The second constraint on our shell is the requirement of 
cubical volumes. The surface area of a shell is roughly given 
by 



S = 4irr 2 = 4n [0.5 • (ri ower + r up 



(10) 



As each shell contains 12 • n% particles, the edge size of one 
piece on the surface can be written as 



^particle — 



12 -n\ 



TV 1 

rT7 (^lower ~~\~ Supper) • (H) 
12 712 
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This edge size should be equal to the width of the shell. 
Setting it equal to d s he\i and solving for 712 leads to 

7f Slower H~~ 



/ /l ' lower \ 1 upper /-i 

ri2 = \— . (12) 

V ±Z ^upper Slower 

While increasing the upper radius of the shell, the second in- 
dex decreases. Since m = and ri2 — 00 for r upP er = Hower 
it is always possible to increase the upper radius until ni 
equals ri2. Having found the upper radius, we place a shell 
of particles at r = 0.5 • (r upP er + Hower) using the coordi- 
nates from the healpix library and continue with the next 
shell. Internal energy and composition of the particle are 
chosen according to the radial coordinate of the particles. 
The initial velocity of all particles is set to zero. 



5.2 Relaxation 

To damp out spurious numerical noise introduced by the 
setup, we relax the object before we start the actual sim- 
ulation. To this end we evolve it for several dynamical 
timescales while applying a time-dependent damping force 
similar to that used by |Rosswog et aL (2004). The specific 
force is given by 



v 2 



(13) 



and applied together with gravitational and hydrodynamical 
forces. The damping timescale r controls the strength of the 
damping term. The smaller it is, the stronger is the damp- 
ing. We start with a large damping force that decreases with 
time and is eventually switched off completely. Afterwards 
we continue the simulation for a short time. If the relaxation 
has been successful, the object will remain in equilibrium 
and particle motions will stay close to zero. Otherwise, as- 
suming that the initial model is not intrinsically unstable, 
the relaxation parameters have to be adjusted. We vary the 
damping timescale with time as 

1 



1 

T 



t < 0.2 t n 



TO 



3(t-0.2 t ma: 
10 °' 6 'max 



0.2 t n 



< t < 0.8 t n 



(14) 



TO 







t > 0.8 t n 



We run the relaxation for a total time t ma x, which is cho- 
sen to be at least several dynamical timescales. The initial 
damping timescale To should be much smaller than the dy- 
namical timescale. Typical choices for a white dwarf mass 
of l.OM are To = 0.002s and t max = 100s. Figure [2] shows 
the configuration of the initial setup and after the relaxation 
step. The final relaxed configuration agrees nearly perfectly 
with the initial profile. Only for the very outermost parti- 
cles of the star the density is overestimated systematically, 
because these particles only have neighbors which are at 
smaller radii and have a higher density. 

5.3 Binary initial conditions 

Having obtained equilibrium models for individual stars we 
need to join two of them to form a binary system. For this 
we use two different approaches, which we will label approxi- 
mate and exact binary initial conditions, following [DanetaT] 



( 2011| ). For approximate initial conditions both stars are set 
on a circular orbit with an orbital period to render the bi- 
nary system marginally unstable. The initial orbital period 
is not determined exactly and can be chosen in different 
ways. Those include an iterative procedure until the system 
is stable for a certain number of orbits or using an approx- 
imation for the distance at the onset of Roche lobe mass 
transfer ( |Eggleton||1983| ). To construct exact binary initial 
conditions we implement the approach used by |Dan et aL] 
( |2011| ). We first set both stars on a stable, circular, synchro- 
nized orbit with a large orbital period such that there is 
no mass transfer and the binary system is stable. Then we 
evolve the binary system in the co-rotating frame and add 
an artificial damping force to remove all residual velocities. 
This is done by adding an additional acceleration a ex t,2 on 
each particle that is given by 



<3-ext,i 



- CJ X (U) X Ti) — 2 UJ X Vi 



Tiamp 



(15) 



where r» and v» are position and velocity of particle i in 
the co-rotating frame, u the orbital frequency of the binary 
system, and Tdamp the damping timescale. The orbital fre- 
quency is updated in each global timestep, i.e. when all par- 
ticles are active at the same time. It is calculated as the av- 
erage of the orbital frequencies necessary to exactly balance 
the total force on the two individual stars. The difference be- 
tween the two individual orbital frequencies calculated for 
both stars is typically smaller than one percent. For the 
white dwarf merger simulations presented in this work we 
set Tdamp to 0.005 s, about a factor of five larger than the 
smallest timesteps in the high resolution run with 1.6 x 10 6 
particles. 

We then apply an artificial radial drift velocity to both 
stars such that the center of mass of the binary system does 
not move. The drift velocity is estimated as described in 



Equation (8) of Dan et al. (2011). We continue the simula- 



tion until the first particle of the less massive white dwarf 
crosses the inner Lagrange point. At this point we switch 
off damping and radial drift, and use the current orbital fre- 
quency to transform the binary system into a non-rotating 
frame. The resulting binary system then gives us the ini- 
tial conditions for simulations of merging white dwarfs as 
presented in Section [6] 



6 APPLICATION TO DOUBLE WHITE 
DWARFS MERGERS 

6.1 Setup 

To study the influence of resolution and initial conditions 
on the properties of the merger simulation we pick a specific 
binary system of a 1.1 Mq primary and a 0.9 Mq secondary 
carbon-oxygen white dwarf. This binary system is of particu- 
lar interest, because it has been shown to lead to observables 
that closely resemble normal SNe la under the assumption 



that a detonation forms during the merger (Pakmor et al 



2012) 



Here we perform four merger simulations of this binary 
system with low (1.8 x 10 5 SPH particles) and high (1.8 x 10 6 
SPH particles) resolution and approximate and exact bi- 
nary initial conditions for each, respectively. The properties 
of these models are summarized in Table [T] The individual 
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Table 1. Model parameters and properties of the four merger simulations. Tabulated are the total number of particles in the simulation 
N p , the type of binary initial conditions used (IC Type), the orbital period T, orbital separation a and total angular momentum L as 
properties of the initial binary system. In addition, several properties of the binary system are shown at the time t when the separation 
is 6 x 10 8 cm. These are the number of particles N^ot on the surface of the primary with a temperature larger than 2 x 10 9 K, the central 
densities of the primary and the secondary white dwarf (p S ec,max and p pr im,max> respectively) and an estimate of the smoothing length 
in the hotspot on the surface of the primary, if present. 
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Figure 3. Initial conditions for the merger simulations. Shown is a density slice in the plane of rotation with a size of 4 x 10 9 cm. 



white dwarfs in these models are first constructed as one- 
dimensional models in hydrostatic equilibrium with a uni- 
form temperature of 10 5 K and composition of equal parts 
by mass of 12 C and 16 C. We then create three-dimensional 
particle distributions resembling the one-dimensional pro- 
files as described in Sec. |5.1| The mass of the particles is 
chosen such that all particles of a particular binary system 
have the same mass. We relax the individual white dwarfs 
for 10 s as discussed in Sec. 15.21 Two white dwarfs are then 
combined to a binary system as described in Sec. |5.3| To 
obtain exact initial conditions we use a relaxation timescale 
of 5 x 10 -3 s and a radial velocity of 1.0 x 10 7 cm s _1 for the 
0.9 M white dwarf and 0.82 x 10 7 cms _1 for the 1.1 M 



white dwarf. We start with an initial orbital period of 50 s 
and stop when the first particle of the secondary white dwarf 
crosses the inner Lagrange-point of the binary system. The 
orbital separation and the orbital period at this time are 
given m Table [l] We then remove the radial and add the 
orbital velocity and restart the simulation in a fixed frame. 

The density distribution of the initial conditions is 
shown in Fig. [3] In the models with exact initial conditions, 
the secondary white dwarf is slightly deformed and fills its 
Roche-lobe. The initial separation of the models with ex- 
act initial conditions is a bit larger compared to the models 
with approximate initial conditions, which also gives them 
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N= 1.8 x 10 5 ,approx. ICs 




Figure 5. Density slice through the binary system when the separation distance is 6 x 10 8 cm. The individual plots are centered on the 
center of mass of the binary system and rotated to have the centers of mass of both white dwarfs on the x-axis. The box has a height 
and width of 4 x 10 9 cm. 



■ N = 1.8 x 10 5 ,approx. ICs - N = 1.8 x 10 6 ,approx. ICs 

■ N = 1.8 x 10 5 , exact ICs - N = 1.8 x 10 6 , exact ICs 




Figure 4. Evolution of the orbital separation for the merger mod- 
els until the binary system has merged. The dashed black line 
shows the separation at which we compare the different simula- 
tions in detail. 



a slightly larger (by ~ 5%) total angular momentum. The 
resolution has no significant effect on the initial conditions. 

The configuration of the Gadget code is mostly iden- 
tical for all four simulations. We use the nuclear reaction 
network described in Sec. [4] for the 13 a-element isotopes 
from 4 He to 56 Ni. It is active only for particles with a tem- 
perature larger than 5 x 10 8 K, because carbon burning is 
negligible at lower temperatures. In the high resolution sim- 
ulations we limit the smoothing lengths of the particles to be 
smaller or equal to 2 x 10 8 cm. This only affects less than 1% 
of the particles, in particular the particles which are ejected 
from the binary system during the inspiral phase. 

6.2 Inspiral 

The main difference between the simulations in the inspiral 
phase is in the stability of the binary systems and thus the 
time it takes until the binary system merges. Fig. [4] shows 
the evolution of the orbital separation for all four simula- 
tions. Obviously, there is a much larger difference between 
models with different binary initial conditions than between 
models with different resolution. In particular, binary sys- 
tems with exact initial conditions are significantly more sta- 
ble than binary systems with approximate initial conditions, 
in agreement with previous results by Dan et al. (2011). In 



our case, the binary systems with approximate initial con- 
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Figure 6. Temperature slice through the binary system when the separation distance is 6 x 10 8 cm. The individual plots are centered 
on the center of mass of the binary system and rotated to have the centers of mass of both white dwarfs on the x-axis. The box has a 
height and width of 4 x 10 9 cm. 



ditions are only stable for 2 to 3 orbits before they merge. 
Their stability shows basically no dependence on resolution. 
In contrast, the binary systems with exact initial conditions 
are stable for about 15 orbits (the lower resolution simula- 
tion) and nearly 20 orbits (the high resolution simulation). 
It is important to note that although the stability of the bi- 
nary systems is quite different the final merger occurs on the 
timescale of about one orbit in all simulations. Note also that 
there are systematic problems with SPH, which can affect 
the mass accretion rate and the stability of the binary sys- 
tem. As mentioned before, the density of the very outermost 
particles of the two objects is overestimated systematically. 
In addition, there is a considerable artificial surface tension. 
Both effects tend to reduce mass transfer and therefore can 
stabilize the binary system. In the end, only high-resolution 
simulations with different hydro dynamical schemes will be 
able to show how long the mass transfer phase before the 
actual merger really takes. 



6.3 The last binary orbit 

The most important stage for the fate of the binary system 
is the final orbit of the merger when the secondary white 
dwarf is destroyed. At this time, it is decided whether a 
detonation forms which leads to the violent merger scenario 



( |Pakmor et al.|2010[ [2011] |2Q12| ) or whether the system ulti- 
mately forms a cold primary white dwarf surrounded by an 
accretion disk and a hot envelope made from the material 



of the disrupted secondary white dwarf (see, e.g. Dan et al. 

[2oTT] ). 

To investigate the effect on resolution and binary initial 
conditions on this phase, we compare all simulations at the 
time when their separation first becomes as small as 6 x 
10 8 cm. This distance is approximately the distance at which 



the system described in |Pakmor et al. (2012) forms its first 
hotspot. 

Fig. [5] shows a density slice at this time. Qualitatively, 
all four simulations are very similar. In all cases, the primary 
white dwarf remains mostly unaffected and is surrounded 
by a very thin envelope. There are, however, some quantita- 
tive differences in the central density of the secondary white 
dwarf which is about to be disrupted. In particular, for the 
low resolution it is about 7 x 10 6 gcm -3 , but only about 
4 x 10 6 g cm -3 for the high resolution simulations. Although 
this difference is dynamically not important, it will change 
the burning products of the secondary white dwarf slightly 
in the violent merger scenario. 

At the place where the material from the secondary 
impacts the primary, a denser region forms. The most in- 
teresting region is located between this dense region and 
the surface of the primary white dwarf. Here, material is 
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compressed and heats up enough to start carbon burning 
which can be interpreted as an indication that a detonation 
can form there. A temperature slice of the binary system 
is shown in Fig. [6] clearly demonstrating the hot region on 
the surface of the primary white dwarf. The plot already 
shows that the temperature is significantly higher for the 
high resolution simulations than for the low resolution runs. 
To quantify this difference, we count the number of particles 
with a temperature larger than 2 x 10 9 K in all simulations. 
We find that the runs with high resolution have the same 
number of "hot" particles within a factor of two (81 par- 
ticles for the high resolution run with approximate initial 
conditions and 44 for the run with high resolution and ex- 
act initial conditions) . The low resolution run with approxi- 
mate initial conditions contains only 8 of those "hot" parti- 
cles and the low resolution run with exact initial conditions 
none at all. Note, that with these small numbers of "hot" 
particles none of the simulations are numerically converged. 



As shown in Pakmor et al. (2011), even with the resolution 



used in the high resolutions presented here the conditions 
in the region where the "hot" particles emerge are actually 
underestimated compared to even higher resolution. 

This result is consistent with other previous simulations 
with low resolution and exact initial conditions that did not 
show any carbon burning particles ( Dan et al.||201l| 2012) 
and simulations with both types initial conditions and high 
resolutions which found "hot" particles ( |Pakmor et al.|2 010 



2011 



Raskin et aT]|2012| ). Therefore, we conclude that the 



resolution of the simulation is the deciding factor which de- 
termines whether hotspots form with particles which started 
carbon burning. The initial conditions, in contrast, do not 
seem to have a significant effect on this phase of the merger. 
This result can be understood in terms of timescales. The 
sound crossing time for the secondary white dwarf with a 
mass of 0.9 M is only about 7 s. Therefore, even with ap- 
proximate initial conditions for which the binary system be- 
comes unstable after a bit more than 2 full orbits, the sec- 
ondary white dwarf has enough time to adjust its structure 
to the gravitational potential of the primary white dwarf. 
Moreover, the difference in the total angular momentum of 
the systems with approximate and exact binary initial con- 
ditions is too small to significantly affect this phase of the 
merger. 



7 SUMMARY AND CONCLUSION 

In this paper we presented the modifications to the Gad- 
get code that are necessary to apply it to problems of 
stellar physics. These modifications include a general equa- 
tion of state, a new efficient nuclear reaction network and 
a mechanism to avoid errors from integrating particles on 
local timesteps. 

We then described a new method to arrange the 
SPH particles in the initial conditions according to a one- 
dimensional density and pressure profile as commonly given 
in problems of stellar astrophysics. We showed that the con- 
figuration we obtain is perfectly stable and retains its radial 
profiles. We then discussed two different ways of setting up a 
marginally stable binary system from two white dwarfs using 



Gadget" code, we investigated the effects of binary initial 
conditions and resolution on the properties of the merger 
of a 1.1 M© and with a 0.9 M© carbon-oxygen white dwarf 
(Pak mor et al.|2012 ). We found that exact initial conditions 
lead to a significantly more stable binary system, confirm- 



ing previous results by Dan et al. (2011). However, we also 



worked out that the type of binary initial conditions has no 
effect on the properties of the actual merger. 

Comparing the simulations in the last binary orbit, 
when they reach the same separation between primary and 
secondary white dwarf, we found that resolution is the dom- 
inating factor. It decides whether carbon burning starts in 
hotspots or not, regardless of the type of binary initial con- 
ditions. We stress that merger simulations should be carried 
out with at least 10 6 particles to avoid incorrect conclusions 
due to under-resolved hotspots. 

Unfortunately, our simulations also confirmed again 
that even for our best-resolved runs the resolution around 
and at the hotspots is far from sufficient to really confirm or 
rule out the formation of a detonation. Given our results it 
seems rather unlikely that this will be possible in the near 
future using an SPH code. Thus, the best way to investigate 
this might be to try to model the accretion stream on a high 
resolution grid as done for helium white dwarf secondaries 



Guillochon et al. (2010). 



approximate and exact ( Dan et al.|2011 ) initial conditions. 
As an example demonstrating the use of the "Stellar 
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